Conformal diagrams for the gravitational collapse of a spherical dust cloud 
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We present an algorithm for the construction of conformal coordinates in the interior of a spher- 
ically symmetric, collapsing matter cloud in general relativity. This algorithm is based on the 
numerical integration of the radial null geodesies and a local analysis of their behavior close to the 
singularity. As an application, we consider a collapsing spherical dust cloud, generate the corre- 
sponding conformal diagram and analyze the structure of the resulting singularity. A new bound 
on the initial data which guarantees that the singularity is visible from future null infinity is also 
obtained. 
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I. INTRODUCTION 



Once a sufficiently large mass is concentrated in a small region of an asymptotically flat spacetime, as occurs in 
the complete gravitational collapse of a star, for instance, a trapped surface appears [l[ and the singularity theorems 
predict that a spacetime singularity forms (see @,| and references therein). An interesting question is whether this 
singularity is naked, or if it is hidden inside a black hole, such that no information escaping from it - light rays in 
i— 1 \ particular - can be detected at future null infinity. 

This question is essentially the contents of the weak cosmic censorship conjecture Q which states that under suitable 
I ' assumptions on the matter fields, the maximal Cauchy development of asymptotically flat, nonsingular initial data 
£~ i . generically yields an asymptotically flat spacetime with a complete future null infinity. In case a naked singularity 
forms, there is a Cauchy horizon, and if it extends all the way to future null infinity it is not possible to predict the 
evolution of a test field in the asymptotic region for arbitrarily large times. As a consequence, future null infinity is 
' incomplete and if this situation persists for perturbations of the initial data, weak cosmic censorship is violated. So 
^ \ far, despite much work, no general proof or disproof of this conjecture has been given. Reviews on this topic can be 
^ '. found in Refs. 

| In this paper we analyze the gravitational collapse in the simple case of the Tolman-Bondi models, describing 
^vq a collapsing spherically symmetric star with zero pressure. In these models spacetime is known in closed, explicit 
form, which greatly simplifies the analysis of its causal structure. In fact, it has been known for a long time that 
such collapse models lead to the formation of shell-focusing singularities, a portion of which is null and visible to 
local observers 0413 • For appropriate initial data, part of the null singularity is even visible to observers which are 
arbitrarily far away from the dust cloud and there is a Cauchy horizon which extends all the way to future null 
infinity. We are interested in understanding how generic this feature is, at least within the class of Tolman-Bondi 
models. 

To this purpose, we develop a method that generates a conformal diagram inside the collapsing, spherical dust 
cloud from given data for the initial density and velocity distributions. This method provides a valuable tool for 
understanding the causal structure of the spacetime. In particular, it enables one to determine in a systematic way 
whether or not a given initial data set results in a singularity that is hidden inside a black hole. Our method is based 
on a combination of analytic and numerical techniques. Analytic tools arc used to understand the behavior of the 
null geodesies in the vicinity of the singularity, while numerical techniques are used to integrate the light rays away 
from the singularities. By generating the conformal diagrams for different initial data sets, we find that it is possible 
to obtain spacetimes with naked singularities which are globally visible without fine-tuning, indicating that these are 
generic within the class of spherically symmetric dust collapse. We also identify a large new class of initial data which 
leads to the formation of such globally naked singularities. 

This paper is organized as follows. In section UH we briefly review the Tolman-Bondi model in the bounded case, 
which describes the complete gravitational collapse of a spherical dust cloud, and state our assumptions on the initial 
data. In section IIII1 we start with a qualitative analysis of the in- and outgoing radial light rays emanating from 
or terminating in the shell-focusing singularity. In particular, we analyze the existence, uniqueness and asymptotic 
properties of such rays in a vicinity of the singularity. Our presentation is self-contained and presents a simple 
derivation of many known results. However, it also goes beyond previous results in the literature insofar that we 
obtain new asymptotic expansions for the light rays terminating at the singularity and a new bound on the initial 
data which guarantees that the resulting spacetime contains a globally visible singularity. Next, in section HVl we 
describe our method for constructing the conformal coordinates inside the collapsing cloud. These coordinates provide 
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a natural extension to the inside of the cloud of the Penrose-Kruskal coordinates for the Schwarzschild metric. Then, 
in section [Vj we present the conformal diagrams corresponding to different initial data, and analyze in which cases the 
resulting singularity is naked or covered by an event horizon. Conclusions are drawn in section IVTl and more technical 
points are discussed in the appendices. 



II. TOLMAN-BONDI DUST COLLAPSE 

In terms of co-moving, synchronous coordinates [13| . the spacetime metric g, four- velocity u and density p for the 
solutions of the Einstcin-Euler equations describing the gravitational collapse of a spherically symmetric dust cloud 
are given by 

g = -dr 2 + ^'51 dR 2 + r(r, R) 2 (<® 2 + sin 2 * d v 2 ), (1) 

where here the function r h-> r(r, i?) describes the evolution of the areal radius along the dust shell R as a function 
of proper time, and r and r' denote the partial derivatives of r with respect to r and R, respectively. We choose R 
such that each dust shell is labeled by its initial areal radius at r = 0, that is, r(0, R) = R. The time evolution of r 
is governed by the one-dimensional mechanical system 

lr(r,R) 2 + V(r(r,R),R)=E(R), V(r,R) := (3) 
2 r 

for each shell i?, where m(R) is the Misner- Sharp mass function [l4j which is determined by the initial density profile 
po according to 

R 



m(R) = 4nG J p {R)R 2 dR, 
o 



with Newton's constant G. The initial data consists of the initial velocity and density profiles vq(R) ■= r(0,R) and 
po(R), respectively, which fix the energy E(R) = v (R) 2 /2 — m(R)/R for each shell R. 

We consider collapsing clouds of finite radius R\ > 0. More precisely, our assumptions on the initial data are the 
following: 

(i) po, v o '■ (0, oo) — > K posses even and odd C°°-extensions, respectively, on the real axis K (regular, smooth initial 
data), 

(ii) po(R) > for < R < R\ and po(R) — for R > Ri (finite, positive density cloud), 

(iii) p'o(R) < for all R > (monotonically decreasing density), 

(iv) 2m(R)/R < 1 for all R > (absence of trapped surfaces on the initial slice). 

Notice that condition (iv) automatically implies that 1 + 2E(R) > for all R > 0, such that equation (fTJ) does not 
exhibit any coordinate singularities as long as r > and r' > 0. Next, we impose conditions on the initial velocity 
profile: 

(v) v (R)/R < for all R > (collapsing cloud), 

(vi) (v {R)/R) 2 < 2m(R)/R 3 for all R > (bounded collapse). 

The condition (vi) means that initially, the potential energy dominates the kinetic one such that the total energy is 
negative, E(R)/R 2 < for all R > 0. One could also consider initial data for which the initial velocity is zero at 
some points, on an interval or everywhere, in which case the data is time-symmetric. In this respect condition (v) 
does not represent a genuine restriction since after an arbitrarily small time such data will evolve into a configuration 
where (v) and (vi) are both satisfied. Nevertheless, our results also apply to the time-symmetric case by substituting 
v a = in our formulae below. Notice that the Oppenheimer-Snyder collapse [l5| . for which E = and the density is 
homogeneous, is not covered by our assumptions. However, this case could be recovered by approximation from data 
satisfying our conditions. 
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For the following, it is convenient to introduce the functions 



c(R) := ,(*) := VE(R)/V(R,R) = J 1 - 

The first quantity is proportional to the the mean density within the dust shell R while q(R) 2 is the ratio between 
the total and initial potential energy. According to assumption (i), these functions have even C°°-extensions on the 
real axis, and assumptions (ii) and (hi) guarantee that c(i?) > 0, c'(i?) < while assumptions (v) and (vi) imply that 
< q(R) < 1 for all R > 0. In fact, our conclusions hold equally well if we replace assumption (iii) by the weaker 
condition: 

(iii)' c'(R) < for all R > (monotonically decreasing mean density). 
Finally, we impose the following two restrictions on the function q. First, 

(vii) q'(R) > for all R > (exclusion of shell-crossing singularities), 

which, together with condition (hi)', implies that r'(R) > for all R > and guarantees that no shell-crossing 
singularities form, see Ref. flo| and the remark below the proof of Lemma Q] in the next section. In order to formulate 
the second condition on q, we first note that the functions c! / R and q'/R are bounded near R = and have even 
C°°-extensions on the real axis. Then, the second condition is 

(viii) For all R > 0, we have q'(R)/R > whenever c'(R)/R = (non-degeneracy condition). 

In particular, this means that the central values of q" and c" cannot be both zero. As we will see in the next section, 
this condition implies the existence of light rays escaping from the central singularity and making it visible, at least 
to local observers. On the other hand, if the central values of q" and c" are both zero, one can show [l(| that no 
such light rays exist. Therefore, the condition (viii) for R = is the key property that determines whether or not the 
central singularity is locally visible. 

Under the assumptions (i)-(viii), the solution of equation © is given by the explicit formula 

r(r,R) = [r 1 (f(q(R)) + V^(RMR) 3 t)}\ (4) 

with the strictly decreasing function 



/ : [0, 1] [0, 7r/2], x n- xyl — x 2 + arccos(x), 



whose derivative is f'(x) = —2x 2 / — x 2 , < x < 1. / is a C°°-function on the interval [0,1). The solution is regular 
on the domain R > and < r < t s (R), where the boundary r = t s (R) describes the shell- focusing singularity which 
is defined by the vanishing of r / R. From equation @, we obtain 

r s {R)= * 2 ZA q{R) \ R>0. (5) 
V4R)q(R) 3 

Since the density p of the dust diverges, the Einstein field equations imply that the Ricci scalar diverges at the shell- 
focusing singularity, and therefore, the boundary points r = t s (R) represent a curvature singularity. The tidal forces 
are much stronger near such points than in the case of shell-crossing singularities Outside the cloud, R > R\, 
the spacctime is isometric to a subset of the Schwarzschild-Kruskal manifold according to Birkhoff 's theorem, see for 



example Ref. [17 1. 



It is worth noticing that the key equation ([3]), describing the dynamics of the dust shells, is identical to its Newtonian 
counterpart if r is identified with Newtonian (absolute) time. What makes the relativistic part much more interesting, 
however, is the analysis of the resulting causal structure of spacctime. A particular interesting question is whether or 
not there exist light rays emanating from the shell-focusing singularity which are able to escape to future null infinity. 
This is discussed next. 



III. LIGHT RAYS EMANATING FROM THE SINGULARITY: QUALITATIVE ANALYSIS 

In this section, we analyze the behavior of the radial null geodesies in the vicinity of the shell-focusing singularity. 
While most of the results derived here are known in the literature, see for instance, Refs. [ij-fllj]. and have been 
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generalized to non-radial causal geodesies [18|, see also [19|, l2C§ for the marginally bound case, our derivation offers an 
alternative, self-contained and simple presentation of the theory. In addition, we obtain new asymptotic expansions 
for the light rays which will be important for generating the conformal diagram in the next section, and new results 
concerning the global behavior of the Cauchy horizon. 

We divide the singular points into the central singularity, So := {( T s(0),0)}, and the remaining part S := 
{(t s (-R), R) : < R < We first establish for each point p £ S the existence of a unique pair of in- and 

outgoing radial light rays terminating at p. Next, we prove that under our assumptions there is a unique ingoing light 
ray terminating at So, whereas there arc infinitely many outgoing light rays emanating from the central singularity 
So. Then, we discuss the global behavior of the Cauchy horizon and determine fairly general conditions on the initial 
data which guarantee that it lies outside the black hole region. In particular, we provide a new upper bound for the 
central density which implies that in the spacetimes developing from initial data satisfying this bound, the central 
singularity is visible from future null infinity. 

Our method for analyzing the light rays is based on new local coordinates (y, i?), where y is defined by r(r, R) = Ry 2 , 
that is, for each point (t, R), y 2 is the ratio between the areal radii of the dust shell R measured at times r and 0. 
With this, the spacetime domain inside the collapsing dust cloud we are interested in is the rectangular region 

A := {(y,R) :0<y<l,0<R< 

whose boundary consists of the union of S U So = {(0,i?) : < R < Ri} (the shell- focusing singularity), := 
{(y, i?i) : < y < 1} (the surface of the cloud), T := {(y, 0) : < y < 1} (the center of the cloud), and II := {(1, R) : 
< R < Ri} (the initial surface), sec figure [1] 
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FIG. 1: A spacetime diagram representing the interior of the collapsing dust cloud in the coordinates (y, R). In these coordinates, 
the interior is the rectangle A which is bounded by the initial surface II, the surface of the cloud O, the singularity EuEo and 
the center V. 



The radial null geodesies are the null curves of the radial part of the metric (p}, 



dr 



r'(r,R) 



dR y/l + 2E{R)' 



(6) 



where e = 1 for outgoing and e = — 1 for ingoing null geodesies. Expressed in terms of the new coordinates (y, R), we 
have 



r _ gm),y) 

and the equation for the radial light rays is 
dy 1 



R 2 



' = y 2 + —y/l-q(R)*y*A(y,R), 

y 



dR 2 



RA(y,R) ( eRQ(R) 



v/1 -<?(#) V) - eQ(R) 



(7) 



(8) 
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where the functions A : [0, 1) x [0, -S- M, Q : [0, -S- M, and g, h : (0, 1) x [0, 1) ->• K are defined as 



c(R) 



y 1 - R 2 q{R) 2 c(R) ' 

f(qy) - /(g) 

9 3 



: = 7= 1 9 /, ^ , , ~ Iff(g'^)- 

Notice that 1 — R 2 q 2 (R)c(R) = 1 + 2E(R) > and c(i?) > are positive, such that Q is well-defined. In the time- 
symmetric case the first term in the expression for A must be dropped, g(q, y) is replaced by the function f(y), and 
the function h which is ill-defined for q — 1 is not needed. It follows from our assumptions that the functions A, Q, g 
and h are C^-differcntiable. As a consequence, the coefficients on the right-hand side of equation ((8j) are smooth for 
all points (y, R) € A in the interior of the cloud. However, equation (|8) is singular at the shell-focusing singularity 
y = and the analysis of light rays emanating from or terminating at those points has to be treated specially. 

Before undertaking this analysis, we summarize the elementary properties of the functions A, Q, g and h which will 
be important later: 

Lemma 1 The junctions A, Q, g and h are strictly positive, C°° -differentiable functions on their domain. Further- 
more, for fixed q, g(q, •) and h(q, •) are strictly decreasing and for fixed R, A(-, R) is strictly decreasing. 

Proof. First, Q > follows directly from its definition and our assumptions. Next, the positivity and monotonicity 
in y of g follow from the fact that / is strictly decreasing. For the function h, we fix < q < 1 and note that 

Oh q 2 y^ 

lim h(q, y) = 0, —(q, y) = . < 0, < y < 1, 

y-*i ay (1 — q z y z ) A l z 

which implies that h(q,-) is strictly decreasing and that h(q,-) > on the interval [0,1). This, together with the 
assumptions (hi)', (vii) and (viii) immediately implies that A > and that A(-,i?) is strictly decreasing for fixed R. 

□ 

In particular, it follows from equation (0 and A > that r' > on A, which precludes the existence of shell-crossing 
singularities, and implies that through each point pG A there passes exactly one pair of in- and outgoing radial null 
geodesies. 

In the following, we analyze the radial light rays through the singular points p € E and p € So- 

A. Asymptotic behavior of the radial null rays near E 

Let p — (0, i?o) G E be a point on the non-central part of the singularity, < i?o < Ri- We show that there is a 
unique pair of in- and outgoing radial null rays terminating at p, which can be parametrized in the form R = Rq + 4>{y) 
for some C°°-function </> : [0,(5) — >• K satisying 0(0) = 0. For this, we first rewrite equation <JSJ> in its inverse form, 



dR 



dy dy V 1 - 1 2 V 2 



RA ( eRQ 



y 2 \ y 



Vi-q 2 y 2 ) -eQ 



-1 



(9) 



In order to find the asymptotic behavior of the function </>, we assume <f> ~ Ay a with i ^ and a > some positive 
exponent to be determined. More precisely, following (Tlj . we set x :— y a and assume 4>{y) = ^{x), with ip : [0, 5) — > R 
a C 1 -function satisfying ^(0) = and dip/dx(0) = A 7^ 0. The equation for the function ip implied by equation (jH) is 



(hp_ _ 2 
dx a 

1 — q^X" 



RA (y - eRQ^/l - q 2 y 2 ^ - eQy :: 



y=x a ,R=R +il!(x) 



Since the left-hand side converges to A 7^ while the expression inside the square parenthesis on the right-hand side 
converges to — ei?oQ(i?o)A(0, Rq) 7^ for x — s- 0, a necessary condition for the existence of such a solution is a — 4. 
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In this case, we obtain, in the limit x — > 0, 



A 2ei?2Q(i? )A(0, J R )' (10) 



In order to give sufficient conditions for the existence of a solution, we make the ansatz 

4>(y)=Ay i [l + z(y)}, 

where A is given by equation (|T0|) and z : [0, 8) — > R is a C°°-function satisfying z(0) — 0. As a consequence of 
equation ([5]), the function z satisfies the differential equation 

dz 

y—+Az = yF(y,z), (11) 
ay 

with the nonlinear term defined as 

yF(y,z) := 



Ay/1 - q 2 y 2 RAy - elVQk^l - q 2 y 2 - eQy* 



According to the definition of the constant A and the elementary properties of the functions q, Q and A, it follows 
that yF : [0,5) x (— Si,6\) — >• R is a well-defined, C°°-function provided that 5 > and Si > are small enough. 
Since yF(0, z) = for all \z\ < Si, it also follows that the function F : [0, 8) x (—81,61) — > R itself is C°°. 

Equation (|11[) is an ordinary differential equation with a regular singular point at y = 0, with the nonlinear 
forcing term yF(y,z). It follows the existence of a unique C°°-solution z : [0,e) — ► K satisfying z(0) = and 
dz/dy(0) = F(0, 0)/5, see Theorem 1 in Q or Theorem [3] in Appendix [Al We summarize the main result of this 
subsection in: 

Proposition 1 Let p = (0,Rq) G S. in t/ie vicinity ofp there exists a unique pair of C 1 radial light rays terminating 
at p. Furthermore, these light rays have the form 

R(y) = R + Ay 4 ll + z(y)}, (12) 

where A is given by equation H10\) , and z : [0, 5) — > R is a C°° -function satisfying z(0) = 0. 

Proof. Existence follows from Theorem[31 as discussed above. As for uniqueness, suppose (y(X), R(X)) is a C 1 radial 
light ray terminating at p. We may choose the parameter A such that the point p corresponds to A = 0. Then, 
according to equation (J9j> , we must have 

dR ,. dR 
lim = lim =Af 0, 

A^o (% 4 ) A^-o 4y 3 d?/ 

which implies that we can choose A = y in a vicinity [0, <5) of p. Furthermore, the C 1 -function z : (0, 8) — >• M, 
^(y) := (R(y) ~ Rv)l(Ay A ) — 1 is bounded and satisfies equation (fTTj) and lim z(y) = according to l'Hopital's rule. 

Now uniqueness follows from Theorem [3] Finally, it follows from equation ([6]) and the sign of A that r increases as y 
decreases to zero, showing that the point p is the endpoint of the light rays. [— j 



B. Asymptotic behavior of the radial null rays near Eo 



Next, we turn our attention to the radial light rays emanating from or terminating at the central singularity, 
p = (0,0) G So- As in the previous subsection, we first try to find the asymptotic behavior by assuming that we 
can write the solution to equation © in the form y(R) = <p(u), where <p : [0,8) — > [0,8i) is a C 1 -function of the 
dimensionless variable u := (R/Ri) a satisfying cp(0) = and dcp/du(0) = A > 0, with a > to be determined. 

The function (p satisfies the differential equation 



dtp _ Ri f— 
du~ 2a V 



q(R) 2 V 2 



RiA(y,R)i 



1 - ei?iQ(i?)u-- 1 - v /l -q(Rfy 2 - eQ(R) 



- y—(p(u),R—Riu a 

(13) 
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The left-hand side converges to A > in the limit u 0. The right-hand side also converges to a positive value if 
a = 2/3, in which case equation (fT3"]) simplifies to 



dl f 3Rl R IU\2 2 



RiMy,R) - 
v 



1 - eR 1 Q{R)y/u-^l - q(R) 2 y 2 ) - eQ(R)Vu~ 

y 



y—ip(u) ,R—R\u^ 



In the limit u — > one obtains 



with 



A := A(0,0) 



Rf^o > 0, 



M %0 , )-^- 5 (g 0) 0)>0, 
g 2c 



(14) 



(15) 



(16) 



where we have defined q := g(0), := q"(0), c := c(0) and c ' := c"(0). 

In order to prove the existence of such light rays, we make the ansatz (p(x) — \x 2 v(x), where x := (R/Ri) 1 ^ 3 and 
v : [0,5) — > [0,Si) is a C°°-function satisfying v(0) = 1. The equation for v follows easily from equation (fl~3|) with 
a = 1/3 and can be cast into the form 



(fa: 



6v J 



9i?i 



TRW 



RiMy,R) 
A 3 



eR 1 Q(R)x 
Aw 



eQ(R)xv 2 



A 



R=R 1 x'5 ,y=Xx 2 v 



Finally, we define a new C°°-ftmction 2 : [0, <5) — > (—(52, 62) through i>(x) 3 = 1 + z(x), such that z(0) = 0. In terms of 
this function we have 



dz 
dx 



6z = xF(x, z), 



(17) 



with the nonlinear term 



xF(x, z) 
9Ri 



-^yi-?(W 



RiHy,R) 



A 3 



eRiQ(R): 



A(l + z)!/3 



^v/i- 9 (i?)V 



s g(i?)x(i + z) 2 / 3 



6. 



_R=K i a:3,y=Aa: 2 (l+z) 1 /3 



Since xF(x, z) vanishes identically for x = according to the definition of A in equation (fTS"j) . it follows that F : 
[0,(5) x (—(52,(52) — * K defines a C°°-function, and we can make use of Theorem |3] to conclude the existence of a 
unique, local solution z : [0,e) — > R of equation (fTT|) such that z(0) = and dzjdxifS) = F(0,0)/7. Since by 
equation ((6]) r increases (decreases) with R for outgoing (ingoing) radial null rays, this demonstrates the existence of 
an outgoing radial light ray emanating from So and the existence of an ingoing radial light ray terminating at So- 
We summarize our findings in 

Proposition 2 Let p = (0, 0) G So &e the central singularity. There exists an outgoing radial light ray emanating at 
p and an ingoing radial light ray terminating at p, which have the form 



R 

Rl 



1/3 



(18) 



y{x) = Xx 2 [l + z(x)}^ 3 , x 

where the constant A is given by equation H15\) and z : [0, 5) — > R is a C°° -function satisfying z(0) =0. 
Furthermore, these rays are the unique C 1 radial light rays y : [0, 6) — > R satisfying lim y(x)/x 2 = A. 

x— >0 

The question of uniqueness is more subtle than in the previous subsection, and deserves a detailed discussion. The 
reason for this relies in the fact that there might exist radial light rays y : [0, S) — > R with y(0) = whose limit 
y[x)jx 2 for x — > does not exist. In order to analyze this, we start with the following technical lemma which is 
proven in Appendix [Bl 

Lemma 2 Let tp : L := (0,6) — > (0,(5i) be a local, C 1 -solution of equation |1^[ ) such that lim (p(u) = 0. Define 

— u— >o 



inf ^ > 0. 



u£l U 



M 



sup < OO . 



Then, the following statements hold: 
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(i) lfm>0, then M < oo. 

(ii) If M < oo and e = — 1, then m > 0. 

(iii) If m > and M < oo, t/ien lim = A, where X is defined by equation U5)) . 

As a consequence of this lemma, we have the following uniqueness result: 

Proposition 3 Let tpa : I := (0,8) — > (0, Si) be the solution of equation of the form tpo(u) = Xuvq(u) with 
vq : I — > M. a C -function such that lim vq(u) = 1, whose existence was proven in Proposition^ Let <p : V := 

u— >0 

(0,5') —> (0, S[) be a local, C 1 -solution of equation |1^| ) with < 5' < S and < S[ < Si such that lim ip(u) = 0. 

— u— >0 
Then, in the ingoing case e = —1, it follows that tp(u) = <fo(u) for all u G that is, tpo is the unique local solution 
which connects Eo. In the outgoing case e = 1, it follows that tp(u) < <fio(u) for all u G that is, ipo is the earliest 
radial light ray escaping from Y,q. 

Proof. Since lim vo(u) = 1 there are constants M > m > such that m < Xvq(u) < M for all u £ I. Let us define 

u— fO 

the function v(u) := Lp(u)/(Xu), u 6 V . According to Proposition [21 it follows that v(u) = Vo{u) for all u £ /' if we 
can show that lim v(u) = 1. 

u— >0 

Now consider the case e = —1 first, and suppose tp(u) 7^ fo( u ) for some u € I'. Since the solutions cannot cross 
on A, it follows that either ip(u) > ipo(u) for all u € I' or ip(u) < <po(u) for all u € I'. In the first case, it follows 
that tp(u)/u > (po(u)/u = Xvq(u) > m > for all u £ I'. Then, Lemma[5]Ji) and (iii) imply that lim ip(u)/u = A, 

u— >0 

and uniqueness follows. In the second case, <p(u)/u < ipo(u)/u = Xvq(u) < M for all u £ I' and Lemma Hfii) and (iii) 

imply that lim ip(u)/u = A. 
u— >o 

Finally, suppose e = 1 and f(u) > tpo(u). Then, Lemma [2{i) and (iii) imply that lim ip(u)/u = A as before, and 

u— ^0 

uniqueness follows. rj 

Notice that in the case of outgoing radial null geodesies emanating from Eo, we cannot prove that <po(u) is unique. 
In fact, it turns out that there are infinitely many outgoing radial light rays emanating from Eo- Indeed, given a 
point p = (R ,y ) S A sufficiently close to E and such that y < (po(u Q ), u = (R /Ri) 2 / 3 , the outgoing light ray 
tp(u) passing through this point cannot cross (po , and so <p{u)/u < <po(u)/u < M for all u G / for which ip is defined. 
Since by Propostion[T]the light ray ip cannot emanate from E, it follows that ip(u) connects Eo- Therefore, there are 
infinitely many radial light rays emanating from Eo, the earliest of which is ipo, which generates a Cauchy horizon. 

We summarize the results obtained so far in the following 

Theorem 1 Given the assumptions (i),(ii),(iii)',(iv)-(viii) made in section 07), we have the following behavior for the 
radial null geodesies near the shell-focusing singularity: for each point p = (0,i?o) G E, i? > 0, there exists a unique 
pair of C°° radial light rays terminating at p. There is also a unique, C°° incoming radial light ray terminating at 
the central singularity (0,0) G Eo. However, there are infinitely many outgoing light rays emanating from Eo- The 
earliest of those is described by the solution given in Proposition^ and generates the Cauchy horizon. 

These results, including the asymptotic expansions (|12I18[) turn out to be important ingredients for the algorithm 
generating the conformal diagrams described in section IIV1 However, before discussing this algorithm, we analyze 
conditions for which the light ray generating the Cauchy horizon arrives at the surface of the dust cloud earlier than 
the event horizon, implying that part of the singularity is visible to outside observers. 

C. Global visibility of the singularity 

Here, we find conditions on the initial data which guarantee the formation of a naked singularity which is visible 
not only to local observers, but also to observers at an arbitrarily large distance from the cloud. This occurs if and 
only if the light ray generating the Cauchy horizon arrives at the surface of the cloud earlier than the event horizon, 
implying the existence of light rays emanating from the singularity which reach the Schwarzschild exterior spacctimc 
at r > 2m. 

In order to analyze this question, it is useful to consider the apparent horizon. In our spherically symmetric 
model, it is determined by the critical surface dividing the two regions in which the areal radius is increased and 
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decreased, respectively, along the outgoing radial null rays. Parametrizing the outgoing radial null rays by R, and 
using equations © and (J7]) we obtain 

* _ *, + ,_ ^ +R H {y ,n } £EmZ^ (l^^E^) (19) 

along the outgoing null rays, where y 2 = r/R, as before. It follows that the apparent horizon is the surface for which 
the expression inside the second parenthesis vanishes, which is equivalent to the condition y = yAH(R) '■= Ry/c(R) or 

r = r AH (R) := R 3 c(R) =2m(R). (20) 

From equation (|T8]). we see that at least for R <C R\ the light ray y = yo{R) generating the Cauchy horizon lies 
outside the apparent horizon, since yo(R) — R 2 ^ 3 while yAn(R) — R- Next, we also note that an outgoing null ray 
emanating from a point inside the apparent horizon, r < 2m(R), cannot escape this region since r decreases while 
m(R) increases as R grows. Therefore, a null ray emanating from a point inside the apparent horizon inside the cloud 
either reaches the singularity inside the cloud, or the surface of the cloud with r < 2m(i?i), and hence it lies inside 
the black hole region. As a consequence, a necessary and sufficient condition for the Cauchy horizon to lie outside the 
event horizon is that the light ray yo lies outside the apparent horizon for all R > 0, that is, 

r {R) := Ry {R) 2 > 2m(R) 

for all R > 0. In the following, we establish that under suitable assumptions on the initial data this estimate holds. 
This is achieved by estimating each term in the propagation equation (|19p for the arcal radius. We start with the 
following simple observation. 

Lemma 3 One has the following upper bound for the function yo, describing the generator of the Cauchy horizon: 

Vo(R) < 1 (j(q(R)) + [f - /(«>)] ^-VE(Rfj =: V (R), (21) 

for all < R < R\, where c(R) := c(R)/c is the normalized mean density profile and where Cq := c(0) and q Q := q(0). 
Proof. Since by equation (|6]) r cannot decrease along yo, we have, using equation (JT)), 



g(q(R),yo(R)) = ^ > ^c(R)t s (0), 

where according to equation ([5]), r s (0) = (tt/2 — f(qo))/{qo\/co)- Now the statement of the lemma follows from the 
definition of the function g and the fact that / is strictly monotonically decreasing. rj 

Notice that rj(R) is strictly positive for R > and satisfies rj(R) ~ i? 2 / 3 for fl<fli. Furthermore, rj only depends 
on the initial function q and the initial normalized mean density c = c/c , but not on its magnitude c . 

Next, assuming that the areal radius is uniformly bounded away from the apparent horizon on some interval, we 
integrate equation (|19l) to obtain an appropriate lower bound on r: 

Lemma 4 Let 8 > ; and suppose yo{R) > (1 + 8)y A H(R) for all < R < Ro- Then, the generator of the Cauchy 
horizon satisfies 

R R 

MR) > + *(l + S)J R 2 <R)dR, m ■■= J R 2 Mv(R), R) ^ ~ ^ >2 ) r,i — dR (22) 



for all0< R< Rq. 

Proof. We estimate the expressions inside the two parenthesis on the right-hand side of equation (fl9|) from be- 
low. For the first expression, we note that y 2 > (1 + S) 2 R 2 c(R), A(y,R) > A(r)(R),R) and y/l - q(R) 2 y 2 /y > 
yjl — q(R) 2 rj(R) 2 /t)(R), where we have used the results from the previous Lemma and Lemma [U For the second 
expression we use 



Pn( m V 71 - ^ pn( p W 1 - {1 + 5 ) 2R2 l( R Y^R) l-(l + 6) 2 R 2 q(R) 2 c(R) 1 1 

RQ(R) y " RQ{R) (ITWW " V 1 TTs - IT 
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Putting everything together, we obtain 

dr 



> 
dR ~ 



(1 + S) 2 R 2 c(R) + R 2 A( V (R), R) V^WviW 



rj(R) 



1 



for all < R < Ro- Integrating both sides of this inequality and observing that r(0) = 0, the statement follows 
immediately. r-j 

Like r/(R), the function is strictly positive for R > and only depends on q and the normalized mean density 
c = c/cq. In contrast to this, the mass function m(R) = R 3 c(R)/2 scales with the magnitude Co of the central density. 
Therefore, under the hypothesis of Lemma 21 we can always arrange for ro(i?) > 2m(R) on the interval < R < Ro 
by rescaling the initial density profile. This observation leads to our final result on the global behavior of the Cauchy 
horizon. 

Theorem 2 Consider an initial density and velocity profile (po,Vo) satisfying 5(0) = 1 and the conditions 
(i),(ii),(iii)',(iv)-(viii) in section\I^ Let p > be sufficiently small such that p < 1 and 

M 2< inf ±4*) (23) 
p " o<fl<i? 1 27 R 3 £(R) v ' 

Then, the reseated data (p Q := p 2 p ,v := pv^) also satisfies the conditions (i),(ii),(iii)',(iv)-(viii) in sectioning and 
the generator of the Cauchy horizon in the resulting spacetime satisfies 



ro(Ri) > ^m(iii) + ^- J R 2 c(R)dR > 2m{R 1 ) 



Proof. First, wc note that £(i?) — R 7 ^ 3 for R -C Ri, and that £(i?) is strictly positive for R > 0, such that the 
function £ i (R)/(R 3 c(R)) is strictly positive and diverges for R —> 0. Therefore, the infimum in the bound (|23|) is 
strictly positive. 

Next, we note that the rescaled solution satisfies c = p 2 c and q = q, since the rescaling is such that both the initial 
kinetic and potential energies scale with p: 2 . Now fix 1 6=1/2 and let Rq £ [0,i?i] be the maximum value for which 
y (R) > (1 + S)y AH (R) for all < R < R . Observe that R > since y (R) ~ R 2 I 3 whereas yAH(R) ^ R for 
R <§; i?i . Now Lemma 0] and the hypothesis imply that 



it 



S f Q 3 /* 

r (R) > -j—jjti 11 ) + S ( 1 + 6 ) R 2 c{R)dR > -R 3 c(R) + - R 2 c{R)dR 

o o 

for all < R < Rq. In particular, this implies that yo(Ro) > 3yAH(Ro)/2. By the maximality of Rq, this means that 
i?o = Ri, and the proposition follows. [— j 

Therefore, given a normalized mean density profile c and an initial ratio q 2 between total and potential energy 
satisfying our assumptions, the existence of a singularity which is locally visible is based on the requirement that 
the second derivatives c and q" at the center R — do not vanish both (sec the comments below assumption (viii) 
in section [TTJ) , while the global visibility is guaranteed if the magnitude cq of the mean density profile c = CoC is 
sufficiently small. For the case of time-symmetric initial data, where q = 1, the expression for the function £ simplifies 
to 

R R 

m = J R lR l{R)2dR = RV{R)2 ~ J r ^ 2dR > = ^ (f V^) 



In the next section, we numerically compute the bound (|23p for a four-parameter family of initial data and show 
that it is consistent with the results obtained from the conformal diagrams. 



1 Other positive values for <5 could also be considered. Our choice is the one that maximizes the constant in the bound 123D . 
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IV. LIGHT RAYS EMANATING FROM THE SINGULARITY: QUANTITATIVE STUDY AND 

CONFORMAL DIAGRAMS 

In this section we present our algorithm for constructing conformal coordinates (T,X), in which the radial part of 
the metric (H|), 



, 2 dR 2 . y/l + 2E(R) 

j{r,Ry r(T,R) 

assumes the simple form 

g = fi(T, X) 2 (-dT 2 + dX 2 ), 

with f2(T, X) > the conformal factor. In these coordinates, the causal structure is transparent since the radial light 
rays are simply described by the straight lines T ± X = const. The coordinates T and X are conveniently obtained 
by introducing the null coordinates U := T — X and V := T + X which satisfy the advection equations 

U = -jU', V = +jV', (24) 

subject to appropriate boundary conditions. Indeed, if U and V satisfy equation (f2"4"| . then 

n 2 (-dT 2 + dx 2 ) = -n 2 dudv = n 2 iiv (-dr 2 + dJ? " 

and the conformal factor is O = 1/ V UV, provided that U and V arc positive. 



A. The special case of a homogeneous cloud with zero initial velocity 

As a simple example, consider first the case of homogeneous density and zero initial velocity, for which q = 1 and 
c = cq = const. In this case, r = Ry 2 , 7 _1 = y 2 / \/l — cqR 2 , where y = / _1 (v / co T ) is independent of R. Therefore, 



y_ 

Co 



2dy 



cadR, 



VI - c Q R 2 



and the substitutions y 
conformal time T, 



cos(T/2), t/cqR = sin(X) lead to the Friedman-Robertson- Walker form of the metric with 



g = 



cos 4 (T/2) , 



CO 



[-dT 2 + dX 2 + sin 2 (X)(d$ 2 + sin 2 d dp 2 )] , < T < tt, < X < arcsin(Vco"i?i). 



The corresponding conformal diagram for Co = 0.75 and R\ = 1 is shown in figure[2] Notice that the apparent horizon 
equation yAu(R) = \fc~oR reduces to the simple equation T = tt — 2X, which describes a time-like three-surface. The 
singularity is spacclikc and hidden inside the black hole region. Surprisingly, the picture changes completely in the 
generic case. As is already clear from Theorem [TJ the central singularity is null and visible to local observers when 
the assumptions (i)-(viii) are met. Furthermore, the null piece of the singularity may extend far enough into the past 
such that the Cauchy horizon lies outside the black hole, see Theorem [2] and the conformal diagrams below. 



B. The generic case 

In the generic case, where the initial velocity and density profiles satisfy the assumptions (i)-(viii), we solve the 
advection equations (f2"4")l by the method of characteristics. The equations (j2"4"|) imply that U and V are constant along 
the out- (t + (R)) and ingoing (r_(i?)) radial null geodesies, respectively, where 
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FIG. 2: Conformal diagram for a homogenous dust cloud with zero initial velocity and the parameter choices cn = 0.75 and 
Ri = 1. The lines denoted by "AH" and "EH" refer to the apparent and event horizons, respectively. Since outside the cloud 
the apparent horizon coincides with the event horizon, the event horizon inside the cloud is determined by the outgoing null 
ray passing through the intersection of the apparent horizon with the surface of the cloud. Notice also that for these parameter 
values, the event horizon intersects the initial surface, implying that the dust particles which are close to the center initially, 
are already imprisoned inside the black hole region. 



In terms of the coordinates (y, R) introduced in section ITTll these equations are 

^±= w+ (y + ,R), ^- =w _ {y _, R) , (26) 

where the function w e (y, R), e = ±1, is given by the right-hand side of equation (|8[). For the numerical integration we 
use, instead of R, a new parameter A and rewrite the equations for the radial null rays as the autonomous systems, 

dy± = ±w±(y±,R±) dR± = ±1 

dX ^l + w±(y±,R±) 2 ' d\ ^l + w±(y±,R±) 2 ' 

The definition of the parameter A is such that it corresponds to the arc length with respect to the (artificial) Euclidean 
metric in the [y, i?)-chart. Therefore, the right-hand sides of equations (|27|) cannot be very large or both very small, 
and there is no need for adaptive or implicit methods when numerically integrating the equations. As a consequence, 
it is sufficient to use a fourth-order Runge-Kutta time integrator (see, for instance, Ref. [2l[) with fixed step size h 
when integrating the ODE system (|27|) . The sign of A is such that r increases with A along the in- and outgoing 
radial light rays. 

Once the characteristics are found, the null coordinates U and V can be constructed by assigning to each outgoing 
null ray a unique value for U and to each ingoing null ray a unique value for V , such that U > 0, V > 0. For this, we 
specify boundary conditions for U and V at the surface O of the cloud, boundary conditions for U at the singularity 
E and boundary conditions for V at the initial surface LI (sec figure [1]). In order to assure that U, V > we require 
that U increases as one moves along O to the future, and then along the singularity S toward the center. Similarly, 
we ask that V increases when one moves along II from the center to the surface, and then along O to the future. 

In what follows, we describe these boundary conditions in more detail. The idea is to match the conformal diagram 
smoothly to the Penrose-Kruskal diagram of the exterior Schwarzschild solution. For a different choice of boundary 
conditions, sec Ref. 1221 ■ 
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Boundary conditions for U and V at 



The quantities U and V at are determined by matching their values to those of the Penrose-Kruskal coordinates 
in the Schwarzschild spacetime. For this, we obtain U and V analytically by computing the trajectory of a free falling, 
radial observer in the Schwarzschild spacetime which starts at r = R\ with initial velocity i>o(i?i). The calculation 
which is performed in Appendix [Cl yields the following result: 



tan(Ue(»)) = — 



1 

Ol 



tan(Ve(y)) 



x exp 

-i( 

oi V 
x exp 





- b\y* - 




\ v 2 


y/i- 




2a\ 


2b\ 





2 A ( 



arctan 




%-bly 2 + yJl-bl 



2a\ 



+ 



2b\ 



l + 2b\ 
bi 



arctan 




where Oi := \/2m{R\)/ R\ = R\^c{Ri) and b± := \J — 2E{R\) = q{Ri)a\. As shown in Appendix [Cl Uq and V@ are 
strictly decreasing functions of y which are normalized such that Uq(0) = Ve(0) = arctan(l) = ir/4. 



2. Boundary conditions for V at IT 



At the initial surface II, we choose V to be a linear function such that the values of it and its derivative match the 
corresponding expressions for the Penrose-Kruskal coordinates. Therefore, 



V U (R) :=V S (1)+P 



R_ 
R~i 



- 1 



P ~Ri 



8Ri 



For the cases studied below we have verified that f3 is positive, implying that V > along II. 



3. Boundary conditions for UafE 

Our choice for U at the spacelike portion S of the singularity is based on the observation that in the Penrose diagram 
for the Schwarzschild spacetime, the curvature singularity at r = is a subset of the straight line 2T = U + V = const. 
Therefore, we choose U such that U + V = const. = Uq{Q) + Ve(0) = n/2. However, this requires the knowledge of 

V at each point of E. We determine the value V p of V at a given point p <= S by computing numerically an ingoing 
radial null geodesies to the past until it intersects either the surface of the cloud <3 or the initial surface IT, where 

V is known. Since V is constant along the ingoing null ray, this determines V at p. Then, U is obtained by setting 
U p := n/2 - V p . 



C. Construction of the conformal diagram 

The boundary conditions determine a unique pair of coordinates (U, V) for each point p G A in the inside of the 
cloud. To find it, it is in principle sufficient to compute a future directed outgoing radial null ray and a past directed 
ingoing radial null ray emanating from p. The outgoing ray intersects either S or 6, where U is given. The ingoing 
ray intersects either IT or O, where V is known. With this algorithm it is in principle possible to construct the images 
of the relevant surfaces IT, T, Q, E, the apparent horizon or a dust particle trajectory in the conformal TX diagram, 
where T = (V + U)/2, X = (V — U)/2, by computing the coordinates (U, V) for each point on these surfaces. 

However, in practice, a problem with this procedure occurs when the outgoing rays approach the singular surface 
E since at E the equation for the null ray is undefined. In order to avoid integrating toward the singularity, we adopt 
the following algorithm: 

Step 1. We choose a uniform grid (yj, i?i) along the surface O of the cloud, where yj = j / J, j = 0, 1, 2, J. Since U 
and V are given at 0, we immediately obtain the conformal image of by computing (£/e(yj), V@(yj)) at all 
gridpoints j = 0, 1, 2, J. 
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Step 2. For each j — J,J— 1, ... wc shoot from (yj, i?i) an outgoing radial light ray to the past. This light ray may 
intersect one of the surfaces II, T, the apparent horizon or a dust particle trajectory at a point p, say. Then, 
we associate to p the value U p = U®{yj) for U, and determine V p by shooting from p an ingoing radial light 
ray directed to the past, as described above. This gives the conformal image of the relevant surfaces in the 
past of the surface O of the cloud. This region necessarily includes a portion of the initial surface II, it might 
also contain all of it plus a portion of the center V, or even all of the center. In the last case, it occurs that 
some past-directed outgoing null ray arrives at the central singularity So for some positive j*, say. According 
to Theorem [T] the subsequent outgoing light rays for j < j* must emanate from Eo as well. 

Step 3. Similarly, we dehne a uniform grid (0,i?fc) along the singularity E, where Rk = Rik/K, k = 1,2, ...,K, and 
shoot from each of these gridpoints an outgoing radial light ray to the past. In order to perform the integration, 
we start at the point (y, R) = (h, Rk + Ah 4 ) close to the singularity, which comes from the truncation of the 
asymptotic expansion (|12jl . Then, we intersect the light ray with the relevant surfaces and determine (U,V) 
at the intersection points like in the previous step. In this way, one obtains the image of the relevant surfaces 
beyond the past of 0. 

Step 4. Finally, wc determine the conformal images of the singularities Eo and E. By construction, the image of E is 
a horizontal line in the TX diagram with endpoints (ir/4,Xo), (7r/4,0), with Xq = (Vo — Uq)/2 defined below. 
The image of Eo, in turn, is defined by a unique value Vo of V and a whole range of values U G [Uo, U±] for 
U. Indeed, according to Theorem Q] there is a unique ingoing radial light ray terminating at Eo while there 
are infinitely many outgoing radial light rays emanating from Eo, the earliest of which generates the Cauchy 
horizon. As a consequence, Eo unfolds into a subset of the straight line T + X = Vo in the conformal TX 
diagram. Therefore, it is sufficient to compute Vo, Uo and U\ in order to determine its endpoints. 

To compute Vo, we integrate the unique ingoing radial light ray to the past, starting from the point (y,R) = 
(1.25ft., -\/4(1.25/i) 3 / (3Ao)) obtained from the truncation of the asymptotic expansion (fT8|) . until it intersects O 
or II. To compute U±, we integrate to the future the outgoing radial light ray generating the Cauchy horizon 
starting from the point (y,R) = (1.25ft, -\/4(1.25ft) 3 / (3Ao)) until it intersects O or E. The intersection point 
determines the value for U±. Finally, the value of Uo is fixed by the requirement that Uo + Vq = 2T = tt/2 since 
it must match the value for U at the left endpoint of E. 

V. NUMERICAL RESULTS 

We apply our method for constructing the conformal diagram to the four-parameter family of initial data 

2 



c(R) := c 



6 / R V 3,„ _ / R x 



'-ri-R-j + 7 (2a - 1} U 



i ^ ) 



where the parameters Co, a, qo and qi are subject to the inequalities < a < 1, < i?jC < 3/2, < qo < 1, and 
< qi < 1 — go- These conditions guarantee that the corresponding initial density and velocity profiles (po,vo) = 
((R 3 c)' j (8irGR 2 ), —Ry/ (1 — q 2 )c) satisfy the assumptions (i)-(viii) on the interval [0,i?i]. Although the resulting 
density profile, 



Po(R) 



8irG 



R\ 2 ^fR^ 



0<R<Ri 



cannot be C°°-smoothly matched to zero at R = R±, it still satisfies po(Ri) = which implies that the metric is twice 
continuously diffcrcntiablc across the surface of the cloud. 



A. Examples of numerically generated conformal diagrams 

Conformal diagrams corresponding to three parameter choices are shown in figures [3j [4] and [5] As anticipated 
in the previous section, the causal structure in the generic case is quite different than in the case of homogeneous, 
time-symmetric collapse illustrated in figure [2] First, we notice from the diagrams that the apparent horizon may be 
spacelikc inside the cloud, whereas it is always timelike in the homogeneous, timc-symmctric case. Next, while Eo UE 
is spacelikc in the latter, E is still spacelikc but the central singularity Eq is null in the generic case, and therefore, 
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FIG. 3: Conformal diagram for the model described in equation (|28|l with the parameter choice Co = 0.1691, a — 0.3, go = 0.75 
and qi = 0.02. The lines denoted by "AH", "EH" and "CH" refer to the apparent, event and Cauchy horizons, respectively. 
In this case, the singularity is hidden inside the black hole region since all light rays emanating from it end at the spacelike 
singularity. Although the center of the cloud appears to be null in the diagram, closer inspection reveals that it is, in fact, 
time-like. The dotted line corresponds to the dust particle trajectory with initial areal radius Ro = 0.95-Ri. The diagram was 
generated with the step size h = 0.0005 and 1500 and 500 points, respectively, for the grids yj and Rk (see steps 2 and 3 in the 
previous section). 

it is visible to local observers. Moreover, the null part of the singularity may either be completely hidden inside the 
black hole region, as in figure [31 or a portion of it may be visible from future null infinity, as in figures 0] and [5] 

It is worth stressing that no fine-tuning is required to construct the examples in which the naked singularity is 
globally visible. In fact, as Theorem [5] shows, it is sufficient to decrease the central density po = 3co/(8ttG) to a small 
enough value in order to produce such singularities. For the initial data described in equation (|28p with the parameter 
values used in figures [3] and 0] the bound ([23]) on Co in Theorem [5] gives 

fi 2 < 0.062, 

which is consistent with the numerical results in those figures. However, the numerical results also indicate that our 
bound is far from optimal, since the transition from local to global visibility occurs around cq ~ 0.169 which is much 
larger than /i 2 . 

We have verified the self-convergence of our numerical results by generating the conformal diagram using different 
step sizes h. In tablc|T]we show the values for the key quantities U\ and Vq in the example of figurc|31 which determine 
the location of the earliest singular point in the collapse, from which the Cauchy horizon emanates. The results 
show self-convergence to an order between three and four. A similar self-convergence test is performed for the dust 
particle trajectory shown in figure[4j corresponding to an initial areal radius of Rq = 0.95i?i. In figurc[6l we show the 
numerical error for different points on this trajectory. The numerical error is estimated by computing the Euclidean 
norm between two successive resolutions (h and h/2, h/2 and h/A and so forth), and clearly decreases as resolution is 
increased. In order to quantify these errors, we show in table [TT] for each fixed resolution the maximum of this error 
with respect to all points on the trajectory. The results exhibit fourth order self-convergence. 
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FIG. 4: Conformal diagram for the model described in equation (|28|) with the same parameter choice as in the previous figure 
except that Co = 0.1688. In this case, there exists light rays emanating from the null part of the singularity which arrive at 
the surface of the cloud earlier than the apparent horizon. Therefore, a portion of the singularity is visible to distant observers 
outside the black hole region. 




FIG. 5: Conformal diagram for the model described in equation (|28|) with the same parameter choice as in figure[3] except that 
go = 0.753. Also in this case a portion of the singularity is visible outside the black hole region. 
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Resolution 


Step size 


Ui 


Error 


CF 


Vo 


Error 


CF 


1 


ft = 0.01 


-5.764905579201557£ 


-01 


— 




— 


1.120836074650528S 


-03 


— 




— 


2 


ft/2 


-5.709420746741624.E 


-01 


5.54848S 


-3 


— 


1.120831798152899£ 


-03 


4.2765.E - 


9 


— 


3 


ft/4 


-5.704917279004125£ 


-01 


4.50347S 


- 4 


12.3205 


1.120831389468409£ 


-03 


4.08684S - 


10 


10.4641 


4 


ft/8 


-5.704566274816609S 


-01 


3.51004B 


- 5 


12.8302 


1.1 2083 1350043764.E 


-03 


3.94246S - 


11 


10.3662 


5 


ft/16 


-5. 704540486 177565£ 


-01 


2.57886S 


-6 


13.6108 


1.120831346181979E 


-03 


3.86179B - 


12 


10.2089 


6 


ft/32 


-5.704538804441139S 


-01 


1.68174E 


- 7 


15.3345 


1.120831345796034.E 


-03 


3.85945B - 


13 


10.0061 



TABLE I: Self-convergence test for the quantities Ui and Vq in the example shown in figure [4] For the i— th resolution the 
error is defined as Ei := \Ui — [7j-i|, and the convergence factor as CFi := Ei-i/Ei. 




FIG. 6: Numerical error for points on the particle trajectory shown in figure [4] The continuous line gives the Euclidean 
distance between corresponding points computed with step sizes ft = 0.01 and ft/2. The dotted one represents the error for 
points computed with step sizes ft/2 and ft/4, and so forth. 



Resolution 


Step size 


Maximum error 


CF 


1 


ft = 0.01 








2 


ft/2 


1.115221660178996£ 


- 7 




3 


ft/4 


6.812854288367279-E 


-9 


16.3693 


4 


ft/8 


4.148846069813869S 


- 10 


16.4211 


5 


ft/16 


2.588085455658281S 


- 11 


16.0306 


6 


ft/32 


1.608422340951133S 


- 12 


16.0908 



TABLE II: Self-convergence test for the dust particle trajectory shown in figure [4] The ith maximum error is defined as 
the largest Euclidean distance between corresponding points on this trajectory computed for resolutions i and i — 1. The 
convergence factor is computed as in table U Notice that a step size of order ft = 0.005 is sufficient to keep the maximum error 
below 1.2 x 10" 7 . 
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B. Phase space diagrams 



As an application of our results, we explore some features of the phase diagram corresponding to the four-parameter 
family of initial data given in equation (|28[) . In figure [7] wc show a subset for fixed q of this diagram, and the critical 
line in the a-co-plane that divides the regions corresponding to initial data giving rise to black holes and globally 
naked singularities, respectively. From the figure, it can be observed that density profiles which are nearly flat close 
to the center, corresponding to a< 1, have a lower critical value for cq than density profiles which are concentrated 
near the center. In this sense, diluted profiles favor the formation of black holes while concentrated profiles favor the 
formation of globally naked singularities. Notice also that although the upper bound for Co obtained in Theorem [5] is 
far from optimal, it nevertheless describes the correct qualitative behavior for the critical line. 

In figure|8]we show the critical line in the a-co-plane for different values of qo and qi = 0.01. As we see, the critical 
value for Co decreases as go decreases from 0.98 to 0.01. Since 1 — q 2 is the ratio between the initial kinetic energy and 
the magnitude of the initial potential energy, this means that large initial velocities in the negative radial direction 
favor the formation of black holes. 




FIG. 7: A cut through the phase diagram corresponding to the four-parameter family of initial data given in equation (|28[1 . 
The cut corresponds to the subset with fixed parameters go = 0.75 and q± = 0.02. The shaded region corresponds to initial data 
giving rise to a naked singularity which is globally visible. The dashed line describes the upper bound for Co from Theorem [5] 



VI. CONCLUSIONS 



We have presented a numerical method for constructing the conformal diagram inside a spherically symmetric dust 
cloud which undergoes complete gravitational collapse. The construction is based on the integration of radial null 
geodesies and a careful analysis of their behavior close to the singularity. Our assumptions on the initial data, namely 
the initial density and velocity profiles, should be physically reasonable and generic within the limitation of spherical 
symmetry and zero pressure. Confirming previously known results [9l-fl2l]. we find that under these assumptions, the 
shell-focusing singularity always consists of a space-like part, which matches the one of the Schwarzschild exterior, 
and a null part which is " naked" in the sense that it is visible, at least to observers which are sufficiently close to the 
center of the cloud. 

Based on our method, we analyzed a four-parameter family of initial data and determined under which circumstances 
the null part of the shell-focusing singularity forms sufficiently early such that a portion of it lies outside the event 
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FIG. 8: The critical line, dividing the regions of initial data giving rise to black holes and globally naked singularities, respec- 
tively, for different values of qo and qi = 0.01. 



horizon, resulting in a naked singularity which is globally visible. We find that this can be arranged fairly simply 
without fine-tuning the data. This shows that it is not difficult to create naked singularities which are in causal contact 
with future null infinity, confirming the results in Q for time-symmetric initial data. More importantly, however, 
these findings also indicate that the formation of a globally naked singularity is stable with respect to perturbations 
within the spherical symmetric, zero pressure model. In fact, this stability statement is confirmed and strengthened 
by Theorem^ which provides a new bound on the initial data guaranteeing that the corresponding spacetime contains 
such a singularity. The corresponding spacetimes are not globally hyperbolic, and as a consequence, it is not possible 
to determine the evolution of test fields obeying hyperbolic partial differential equations with light speed propagation 
beyond the Cauchy horizon without imposing boundary conditions at the null singularity. Since the Cauchy horizon 
is located outside the black hole region, this implies, in particular, that it is not even possible to predict the evolution 
of such test fields in the Schwarzschild exterior for r > 2m. 

Besides the results already mentioned, there are several possible applications of our numerical algorithm. First, it 
provides a tool to systematically analyze the phase space of spherically symmetric dust collapse, identifying the class 
of initial data leading to globally naked singularities. As an example, we have explored in section fV Bl different cuts 
through the phase space, and determined the critical set within these cuts which separates the black holes from the 
globally naked singularities. It should be interesting to perform a more exhaustive analysis, including, for example, 
initial data with arbitrarily concentrated density profiles or families of initial data which contain the Oppenheimer- 
Snyder model as a limiting case. Second, the numerical construction of the conformal coordinates may be useful for 
describing the propagation of test fields on the collapsing spacetimes considered here. For example, the wave equation 
= reduces to a family of flat space wave equations of the simple form 

" ^ + V(T,X)i/> tm = 0, $(T,X,0,<p) = ^J2^rn(T,X)Y im ^, l p) (29) 

with a potential V, when decomposed into spherical harmonics Y tm . Therefore, properties of the solutions outside 
the Cauchy horizon can, in principle, be inferred from the structure of the potential V. 

It remains to be seen whether or not our results can be extrapolated to the realistic gravitational collapse, in which 
pressure and angular momentum may delay the formation of the singularity. If these effects are sufficiently strong, 
it is imaginable that the resulting singularity is hidden inside the black hole, such that weak cosmic censorship is 
valid. Unfortunately, analyzing scenarios with more realistic equations of state in the absence of symmetries requires 
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much more efforts since exact solutions are not available anymore, at least not for generic initial data. Intermediate 
steps toward understanding the general case are the following: first, it should be possible to generalize our method 
to the case of a spherically symmetric collapse with nonzero pressure, since it is based on the numerical integration 
of radial null rays and their asymptotic behavior near the shell- focusing singularity. Next, one could consider a 
collapsing fluid star which is slightly nonspherical, in which case the problem may be analyzed using perturbation 
theory for spherically symmetric spacetimes. In the dust case, numerical work [23| and recent analytic results [Hf for 
the self-similar collapse show that linear fluctuations from smooth initial data grow without bound when approaching 
the Cauchy horizon, providing evidence that the latter is unstable at the linearized level. Finally, the problem can 
be approached using numerical approximation. In [25j . a numerical code was used to evolve collisionless prolate 
gas spheroids, and it was shown that when sufficiently large, these spheroids form a curvature singularity before an 
apparent horizon appears. However, as pointed out in (26|, this does not exclude the formation of trapped surfaces in 
the maximally extended spacetime, and so the numerical work in [25j does not imply a violation of cosmic censorship 
at this stage. More recently, the gravitational collapse has also been analyzed in five [U HH and higher-dimensional 
spacetimes p9ll30l|. 

Even if it turned out that generically, singularities in the nonspherical gravitational collapse with a realistic equation 
of state are hidden inside black holes, it would still be interesting to understand what happens in the limit when the 
collapse becomes spherical and pressure can be neglected. Presumably, the fact that in this limit naked singularities 
which are visible from future null infinity exist should have an imprint on the perturbed case. 
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Appendix A: An existence theorem for nonlinear perturbations of regular singular points 

In this appendix, we include a short proof of the following theorem which we used in section IIIII to show the local 
existence of light rays emanating or terminating at the singularity. It is based on basic tools from the theory of 
dynamical systems. 2 For a generalization to systems of equations and applications to relativistic stars, see Ref. |31| . 

Theorem 3 Let f3 > 0, and let D C R 2 be an open subset o/R 2 which contains the origin. Furthermore, let f : D — > R 
be a C°° -function. Then, the differential equation 

x^-+0y = xf(x,y) (Al) 
ax 

has a unique local C°° -solution y : (0,e) — > R which is bounded. Moreover, this solution satisfies 

v n n ,. dy /(0,0) 
hm y(x) = 0, hm — (x) = — — . 

x^o yy ' ' x-vo dx y 1 /3 + 1 

Proof. Define a := /(0, 0), and introduce the parameter t — — log(x) for x > 0. Then, the solutions of equation (|A1[) 
are given by the trajectories of the autonomous system 

= Au + F{u), (A2) 



with 




U(M) -/(«,,), ' <A3) 



The dynamical system described by equation (|A2[) has a stationary point at u = 0, and the linearized system at 
this point is given by the matrix A. The eigenvalues of A are —1 and /3 > 0, with corresponding one-dimensional 



2 We thank Thomas Zannias for pointing out to us the elegant proof of this theorem. 
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stable and unstable manifolds describing the sets of points in phase space converging to u — for t — s- oo and 
t — > — oo, respectively, see for example Ref. [32|. Therefore, the solutions of equation (| Al|) which are bounded for 
small x = cxp(— t) correspond to the stable manifold. Its tangent vector at u = is given by the eigenspacc of A 
corresponding to the eigenvalue — 1. Since this eigenspace is generated by the vector ((3 + l,a), the slope of the 
solution at u — is a/(/3 + 1). rj 



Appendix B: Proof of Lemma [2] 



Here we prove Lemma [2] which is a statement about local, C 1 -solutions (p of equation (|14[) satisfying lim ip(u) = 0. 

u— >0 



For this, it is convenient to rewrite equation (|14|) in the form 

2 



— = A(y,u) - 



da 



1 - eB(y,u)y/u - 

y 



- ^ B (y, u )Vu, 



(Bl) 



where the functions A, B : [0, 1) x [0, 1] —¥ M arc defined as 

A(y,u) := -RlK{y,R)^l-q{R)^ 



R=R 1 u 3 / 2 



B{y,u) := R!Q(RWl - q(R) 2 y< 



R=R lU 3 / 2 

According to Lemma [TJ these functions are continuous, and they satisfy 

A Q := 4(0,0) = ^i??A(0,0) = A 3 > 0, 
B Q := B(0,0) =i?iQ >0, 

where A is defined in equation (fT5)) and Qo := Q(0) > 0. Therefore, given S > with 6 < min{A , B }, there exists 
6 y > and S u > small enough such that D := [0, S y ] x [0, S u ] C [0, 1) x [0, 1] and 

0<A Q -6< A{y,u) <A Q + S, 
< B a - S < B(y, u) <B Q + S 

for all (y, u) € D. Let ip : (0, 5 U ) — > (0, S y ) be a local solution of equation (|F31[) such that lim <p(u) = 0, and set 

u— >0 



inf ^>0, 



M 



sup < oo. 



0<u<<5„ U 0<u<<5„ U 

We are now ready to prove the lemma. 

(i) Suppose m > 0. Then, we can use u/y < 1/m in equation (|L31[) and estimate 



t , B + 6 nr 

1 H vo„ 



3(B + ^) 



M', 



for < u < 5 U . Since lim y(u) = this implies that ip(u) < M'u for all < u < 8 U . Therefore, M < M' < oo. 



(ii) Conversely, suppose M < oo, and assume e = — 1. Using the estimate u/y > 1/M in equation (|L31[) and the 
positivity of _B on Z3, we obtain 

dy A - S 



> 



=: rri > 0, 



du ~ M 2 

which implies tp(u) > m'u for all < u < S u . Therefore, m > m! > 0. 
(iii) Finally, suppose m > and M < oo. According to the L'Hopital's rule (see, for instance, Ref. HH), we have 



lim ^ } = lim 



u— >o u 6 



u-¥0 \ u 



ip{u)\ dip 3 

— — — A , 



where we have used equation (|F31[) and the boundedness of tt/y in the second step. Therefore, lim ip(u)/u = A. 

u— >0 

□ 
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Appendix C: The surface of the cloud in Kruskal coordinates 



The purpose of this appendix is to construct Penrose-Kruskal coordinates along the surface of the cloud, 9. For 
this, we start with the expressions for the Kruskal null coordinates which are related to the standard Schwarzschild 
coordinates (t, r) by p|: 



U 
V 



± 



r 
2m 



exp 



t 



r 

2m 



exp 



4m 

r + i 
4m 



where the choice of sign (+/— ) corresponds to the region inside or outside the event horizon, respectively. These 
coordinates satisfy the relation 



UV = 11 — 



2m 



exp 



r 

2m 



r > 0. 



(CI) 



The surface is generated by the trajectories of freely falling particles with zero angular momentum in the 
Schwarzschild spacctime, for which the equations of motion are given by equation ([3]) with R = R±. The coordi- 
nate t along the geodesic is determined by the conservation of energy equation 



dt _ yT + 2E 1 
dr ~ 1 - 



2 m i 



where E\ := 2E(R±) and mi 
equations, 



Using this, equation ((3j) and the definition of U and V, we find the following 



4 TOl Alo g ([/) = 1 

dr r — 2rai 



4mi— log(V) 
dr 



1 



1 + 


2E l \ 


2m i 

r 


+ 2E 1 ) 
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2E 1 \ 


2mi 
r 


+ 2E l J 



r — 2mi 

which are valid both outside and inside the event horizon. Integrating equation (|C2[) . we obtain 



(C2) 
(C3) 



U(y) 



1 



if - b\y 2 -yJl-bl) exp 



y 

2a\ 



2b\ 



l + 2b( 
bi 



arctan 



hy 



where a\ := 2m\/R, b\ := —2E\ and y 2 = r/R\. From the relation (|Clj) it follows immediately that 

'l + 2b\ 



V{y) = - 
ai 



&?y 2 + 2/ v / i-^)cxp^ 



" +.v /r ^? 
i 



2b\ 



arctan 



hy 



y 'a 2 



b\y 2 



A translation in t induces the transformations U i— > nil and V i-> V/k, with n > 0. We choose k such that 
U(0) = V(0) = 1. The Penrose-Kruskal coordinates Uq and Ve in section HVl are obtained from the resulting rescaled 
quantities after compactification, f/e(y) : = arctan([/(j/)), Ve(y) := arctan(F(y)). 

Using equations (|C2[) and (|CJ3[) and the fact that r < along the surface, it is not difficult to check that U and V 
are positive everywhere along 0. Therefore, the conformal factor in section HVl SI = 1/ \/U@V@, is well defined. 
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